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Abstract 

In this article a generalized mathematical model describing the interactions between malignant 
tumour and immune system with discrete time delay incorporated into the system is considered. 

Time delay represents the time required to generate an immune response due to the immune 
system activation by cancer cells. The basic mathematical properties of the considered model, 
including the global existence, uniqueness, non-negativity of the solutions, the stability of steady 
sates and the possibility of the existence of the stability switches, are investigated when time delay 
is treated as a bifurcation parameter. The model is validated with the sets of the experimental data 
and additional numerical simulations are performed to illustrate, extend, interpret and discuss the 
analytical results in the context of the tumour progression. 

Keywords: delay differential equations, stability analysis, Hopf bifureation, immune system, 
tumour, immunotherapy 


Acknowledgements 

This work was supported by the Polish Ministry of Seienee and Higher Edueation, within the Inventus 
Plus Grant: ’’Mathematieal modelling of neoplastic processes" grant No. IP2011 041971. 









2 


1 Introduction 


Cancer immunology is a new trend in immunology that foeuses on the variety of interaetions between 
the immune system and eaneer eells. The aim of eaneer immunology is to diseover new eaneer 
immunotherapies to fully eure or at least to retard the progression of the disease. Aeeording to reeent 
reports [[T| lymphocytes sensitive to tumour antigens are often present in the patient’s body, providing 
a potentially natural way to prevent tumour growth. On the other hand, it is also observed that the full 
immune response usually does not occur and so the natural proteetion afforded by the immune system 
is ineffective. 


In general the interaetions between tumour cells and immune system are very eomplex and involve 
a number of different kinds of molecules, immune eells and events, 0- The immune response against 
tumour antigens involves different kinds of immune eells: B-Lymphoeytes producing and secreting 
antibodies into the blood or placing them on their surfaees, eytotoxic T-Lymphoeytes - the effec¬ 
tor cells that destroy the antigens and thus the antigen-bearing eells, T helper-lymphoeytes seereting 
interleukins, hence stimulate both T and B eells to divide and finally T-suppressor cells that end the 
immune response. All these eells are produced in the body from a small fraetion of the preeursor cells. 
Whenever antigens are presented to the preeursor eells they are aetivated, start to proliferate and then 
differentiate into effector or memory cells. Unfortunately, these mechanism do not guarantee the fully 
effieient response of the immune system due mainly to the tumour’s own defenee system. Tumour 
eells proteet themselves from eontaet with lymphoeytes by, amongst other methods, surfaee expres¬ 
sion of ligands whieh initiate the apoptotie signal in the eytotoxie T-Lymphoeytes, thus killing off the 
eytotoxie T-Lymphoeytes before the neeessary eell-eontaet with the tumour eells ean be initiated. 

Clearly, due to the eomplexity of the considered problem systems describing most of the tumour- 
immune system interactions would involve a large number of variables and thus equations, as in Q, 
where a model of 12 equations is proposed and numerieally analysed or in Q, where a system of 5 
equations deseribing the response of the effeetor eells to the growth of tumour eells is studied. More 
reeently, systems of four or more equations investigating the immunotherapy and its improvement 
for high grade gliomas ( [|^[6|), predieting outeomes of prostate eaneer personalized immunother¬ 
apy ( Q) or humoral mediated immune response [[^ were proposed. Also well known immune 
system Marehuk’s model, fflTT), deseribing the immune system dynamies which was applied to 
the description of the tumour growth phenomenon, see e.g. [ 12|. However, to better understand the 
main underlying meehanism of immune response to the presence of the tumour eells eomplex models 


were simplified and thus redueed into two- or three-equation and less detailed ones, see e.g. [[TT 
or 


191 


In this article the influenee of a particular modification to the generalization of the two-equation 
model proposed in p9| is studied. The proposed generalization eoneems the form of the stimulus 
function, namely only some essential assumptions regarding the shape and regularity of this function 
are assumed for the theoretieal part of the paper. The proposed modification reflects the faet that the 
immune system requires some time before the presence of a tumour (aetually, tumour antigens) in the 
body is deteeted and the immune system is able to respond. This assumption leads to the ineorpora- 
tion of a diserete time delay in the funetion deseribing the immune system stimulation by presence of 
the tumour eells. For the obtained system the basie mathematieal properties of the model, the possible 
asymptotie behaviour of the system depending on the delay and the role of the native background im¬ 
munity and the immune stimulation strength are investigated. For the partieular form of the stimulus 
funetion the eonsidered model is validated with two sets of the experimental data for partieular tumour 
eell line and further studied in the eontext of the model dynamies within the framework of the other 
experimental data. It is shown that model ean refleet different patterns, ineluding periodic, reported 
in the experimental literature. For the estimated parameter values it is shown that an increase of the 
baekground immunity or the stimulation strength of the tumour cells on the immune system enlarge 
the stability region for the positive dormant steady state. Moreover, the increase of the baekground 
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immunity might stabilize the unrestrieted tumour growth observed for the eontrol experiment. 

The paper is organized as follows: in Section the model without delay and summary of the 
existence and stability of the non-negative steady states are presented; in Section the model with 
discrete delay is described and analysed. Next, in Section]^ model is validated with the experimental 
data. The numerical simulations confirming and extending the analytical results are shown in Section 
Finally, in Section above results are summarised and discussed. 


2 Presentation of the base model 


In [191 a simple model for tumour-immune system interaction, that takes into account the two most 
important variables: the size of specific anti-tumour immunity at time s (X(s)) and the size of tumour 
(T(j')), is studied. To describe the size of the overall specific immunity against the tumour Forys et 
ah, (Tg, make the following assumptions. First, in the absence of tumour antigens one observes a 
constant (with rate w) production of precursor cells that are able to respond to the tumour antigens 
and the natural death of the cells, which is expressed by a linear term uX. Thus, as a consequence 
of the tumour absence within the body a constant immunity level (so-called background immunity) is 
kept. Second, the presence of tumour antigens stimulate an increase of the immunity level. Moreover, 
it is assumed that immunity response is proportional to the current size of the immunity X with the 
proportionality coefficient (aF), which is assumed to be a bounded function of both immunity size and 
tumour antigens or tumour antigens only. It is also assumed that the creation of the antibody-antigen 
complexes lead to the death of the immune cells and it is modelled by a bilinear term XY with constant 
coefficient b. Clearly, if for particular tumours this process is not observed, then b = 0. Finally, the 
destructive influence of immune system on the tumour development is described by the bilinear term 


XY with the constant rate c, similarly as it was assumed in [ 16|. Additionally, it is assumed that the 
tumour grows according to the exponential law, as e.g. in [14,T^20|, with growth rate equal to r. 
The system proposed in [[Tg has the following form 


dX 

d5 


w-uX + aF(X, Y)X - bXY, 
Yir - cX), 


(2.1a) 

(2.1b) 


where the presence of tumour antigens that stimulate the immune system response (F) is defined by 


Fi{X,Y) 


ytT 

+ 7 “’ 


or F 2 iX,Y) 


k" + 7“ ■ 


( 2 . 2 ) 


Different forms of Fi and F 2 reflect different assumptions on tumour-immunity interactions, and both 
of them were earlier considered in the literature in the context of the immune system response. For 
F = F 2 it is assumed that the number of stimulated cells does not depend on the size of the immune 
system itself, hence stimulation is signal dependent. Such a form of the stimulus function in the 
context of immune response was considered for example in [[4 27 -29| and more recently in [ TgiTj 
25|. In case of F = Fi a number of signal molecules need to be reached before faster production 


of immune cells is activated (after which, the size of anti-tumour immunity increases) hence such a 
function is often called a ratio-dependent stimulation function and it was considered in the immune 
system - tumour interactions context in [|g. 

Due to the biological interpretation of constants, it is assumed that all parameters are non-negative 
and in particular parameters a, c, u, r, w and k,, i = 1,2 are strictly positive. Moreover, it is assumed 
that a, which describes the switching characteristics of the stimulus functions, is greater or equal to 
1 . 
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In order to simplify the model analysis, in the analytieal part of the paper, its non-dimensional 
form is eonsidered 


= d(a> - X + pf(x,y)x - if/xy), (2.3a) 

= y(l-.r), (2.3b) 

where to = cwl(ur) > 0 is the native (baekground) immunity, i.e. the natural level of immunity in the 
absenee of antigens, p = aju > 0 is the stimulation strength, ifj = bY^ju > 0 (Fq = kir/c or Fq = ^2 
depending on the form of F funetion eonsidered) is the strength of anti-immune aetivity, d = u/r > 0 
and f(x, y) is the stimulus funetion given by 


dx 

dt 

dt 


Mx,y) 


X" -l-y“ 


or f 2 {x,y) 


1 +y“'’ 


O' > 1, 


(2.4) 


time and variables are sealed t = rs and x = cXlr,y = F/Fq, respectively. 

As mentioned before in the well established literature the stimulating effect of the presence of 
tumour antigens on the immune system response is usually described by the function f = f 2 and often 
( PPH [^[^ ) with a = I (that is by the Michaelis-Menten function). Hence, that form is considered 
in the part of paper devoted to the model fitting to the experimental data and its further validation. 


However, function /2 given by ( |2.4[ ) has several general properties, thus for the theoretical part on the 
paper a more general family of stimulus functions, fulfilling conditions: 


(Al) the function / is a non-negative Lipschitz continuous on R> = [0, -too); 
(A2) f(y) is bounded by one for all y 6 R>; 

(A3) /(O) = 0, lim fiy) = 1, /(I) = 1/2; 

y —>+00 

(A4) f(y) is increasing for all y 6 R>; 

(A5) /(y) is differentiable for all y 6 R>; 

(A6) the function / has at most one inflection point, 
and if further needed the additional smoothness 


(A7) /(y) 6 for all y 6 R> 


is considered. 

Clearly, assumption (A[^ together with the form of the system ( |2.3[ ) guarantees existence of non¬ 
negative unique solutions of ( |2.3[ ) with non-negative initial condition ( |3.2[ ). Additionally, the global 
existence of the solution is a consequence of the assumption (A[^. 

In [ p[9| the importance of parameters to (background immunity) and p = to -1 +p, (so-called overall 
immune capacity) is stressed since both of them strongly influence the existence and the stability of 
the steady states. Hence, due to the biological interpretation of the parameter p in the rest of the paper 
its positivity, i.e., 

m -I- p > 1 

is assumed. 


Because of the form of Eq. ( |2.3b[ ) all non-negative steady states are always of the form 0) or 
(1,^), where ^ > 0. Assumption /(O) = 0 (see (A[^) implies that there always exists a semi-trivial 
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steady state A = (a>,0). Clearly, the existenee and number of positive steady states depends on the 
number of solutions of equation 

pf(y) = >ffy + l-(o. (2.5) 


Assumption (A|^ implies that there exist at most three positive solutions of ( |2.5| ). If ^ = 0, then 
due to (A[^ and (AQ there exist no positive solutions of ( |2.5[ ) for oj > I and there exists exaetly one 
yR for max{0,1 - p} < m < 1. If i/r > 0, then the situation is more eomplieated. If the funetion / has 
no infleetion point (the assumption (A[^ implies that / is eoneave), then for 0 < cu < 1 Eq. ( |2.5| ) has 
zero or two solutions in a generie ease. On the other hand, for cu > 1 Eq. ( |2.5| ) has exaetly one positive 
solution. If the funetion / has one infleetion point, Eq. ( |2.5| ) has one or three positive solutions for 
OJ > \ and zero or two for 0 < m < 1 (in a generie ease). 

Arguments analogous to those presented in imply that the stability of the steady states of the 
system ( |2.3| ) with / fulfilling eonditions (A[^-(A[^ are the same as for particular forms of / functions 
given by ( |2.4[ ). The results regarding the existence and stability of the steady states of ( |2.3| ) are 
summarized in Tab. [J and partially in Eig. Note that, in additionally it is shown that there are 
no periodic solutions of the system ( |2.3| ). 


Table 1: Existence and stability of the steady states of the system ( |2.3[ ) for / fulfilling conditions 
(A[^-(A|^ depending on values of parameters a, oj, if/, p. Numbers in the column ips indicate the 
number of possible inflection points of the function /. 


if/ OJ ips 

Steady states 

Stability 

Semi-trivial steady states 

if/>0 

0 < 0) < 1 

0,1 

A = (w,0) 

unstable 

OJ > 1 

locally stable 

Positive steady states 

i/r = 0 

o<oj<r 

0,1 

R = {l,yR),yR > 0 

globally stable 

OJ > 1 

none 

if/>0 

0 < OJ < 1 

0,1 

none 

B^(l,yB), C = (l,yc) 

0 <yB <yc 

B locally stable 

C unstable 

OJ > 1 

0,1 

C = (l,yc),yc > 0 

B unstable 

1 

B^(l,yB), C-(l,yc), £> = (hyn) 

0 <yB <yc <yD 

B unstable 

C locally stable 
D unstable 


* under assumption oj + p > \ 


Since in a later section of the paper / is given by the Michaelis-Menten function, the conditions 
for the existence of positive steady states of the system ( |2.3[ ) for a = 1 and f = f2 are directly provided 
below (compare with Eigure 1). 


Proposition 2.1 Let if/ = 0 and f = f 2 is given by (2.4), then 


(i) for OJ > \ there are no positive steady states of the system ( |2.3[ ),' 

(ii) for maxjO, 1 - p} < m < 1 there is exactly one positive steady state of the system (|2.3|) given by 


R 


\ - OJ 


(!,»)= 

\ ^p + OJ - I' 


Proof : Because of the form of Eq. ( |2.3b[ ) all positive steady states are always of the form (1,^), 
where ^ > 0. Clearly, the existence and number of positive steady states depends on the number of 
solutions of p.5| ). Eor if/ = 0, due to the facts: /(O) = 0, lim f(y) = 1 and / being increasing, there 

y—^+oo 

exist no positive solutions of p.5|) for oj > I and there exists exactly one ^ for maxjO, 1 -p] < oj < 1. ■ 
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Proposition 2.2 Let if/ > 0, a = I, f = f 2 is given by ( |2.4[ ) and A = (l + (l-6t; - pflij/Y -A{\-oj)lij/, 
then 


(i) for oj > \ there exits a unique positive steady state C = (1, jc) of the system (2.31 with 


yc 




( 2 . 6 ) 


(ii) for 0 < < 1 and p > { vjij/ + Vl - cof two positive steady states (I.Jb) ond (l,3^c) of the 

Wp + "-i_,_ vaI 


system (2.31 exist where 


ys = 




and yc is given by p.6| ). On the other hand, ifO<co< 1 and p < {v/f + Vl - cu) , then there 
are no positive steady states of the system (Q. 


Proof; 


Eq. ( |2.5| ) takes form py = [fy + 1 - + y). Dividing both sides by if/ and denoting [5 = and 


V = 7 one obtains 


y^+ {\-v+p)y+f = 0. 


(2.7) 


Inequality to > \ implies f3 <0, henee Eq. (2.7) has one positive solution and the part (i) follows. 


Eor 0 < m < 1 there is jS > 0. If real roots of ( |2.7| ) exist, then they have the same sign. Moreover, 
they are positive if and only 

V>l+yS. (2.8) 


On the other hand, the existenee of real solutions to ( |2.7[ ) depends on the sign of the determinant 
A = (1 - V + - 2v(l +/?) + (! -fSf. It ean be easily ealeulated that A > 0 if and only if 


< (1 - or y > (1 + 


(2.9) 


Note, that foryS > 0 the first inequality of ( |2.9[ ) is equivalent to y < 1 +f3-2 s/f, whieh eontradiets ^ 

On the other hand, the seeond inequality of (|2.9|) is equivalent (again for /3 > 0) to y > 1+^ + 2 y/f. 


whieh implies (2.8 1 . Thus, the part (ii) follows. 


Y = 0, a = 1 


Y > 0, a = 1 




Figure 1: The to-p spaee divided into regions where different numbers of steady states for cr = 1 and 
if/ = 0 or iff > 0 are observed. Numbers and text in the regions indieate the number of positive steady 
states of the system ( |2.3| ) for f = f 2 and theirs stability for r 
stand for stable and unstable, respeetively. 


0; abbreviations stab, and unstab. 
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3 Model with discrete time delay 


The immune system requires some time before the presence presence of the tumour (or rather, the 
tumour’s antigens) in the organism is detected and the immune system is able to respond. From the 
mathematical point of view, to reflect such phenomena, one can include time lags into the system, 
e.g. by introducing a discrete delay into the stimulus function. Although the model itself differs 
from the one considered in pT| , the idea presented there of introducing one discrete time delay only 
in the function describing the stimulating effect of the tumour on immune cells proliferation in the 
system variable representing the tumour cells compartment is followed here. That assumption reflects 
situation when the immune system response due to the presence of the tumour cells depends on the 
amount of the immune system cells at a given time t and on the amount of tumour cells at time t - r, 
where r > 0 is the time needed for the immune system to develop a suitable response after antigens 
recognition. Following that approach time lag is introduced in a classical way, as it is done in the 


logistic model by Hutchinson [30| or in the Gompertz model [3F32|, where delay is introduced to 


the function representing the per capita growth rate. The same assumption of the form of time delay 


in the model of the immune system-tumour interactions has been recently made e.g. in [24|. Clearly, 


in the literature 115 17 181 discrete delay in the variable representing the size of the immune system 
was also considered. However, this idea is not followed since the priority is given to the immune 
system’s time reaction. It is also assumed that the immune system does not require time to know the 
size of itself. 

The non-dimensional system with one discrete delay has the following form 


dx{t) 

dt 

dy(0 

dt 


= d{a) - x(t) + pfiyit - T))x{t) - il/x{t)y{t)), 
= y(t)(l -x{t)). 


(3.1a) 

(3.1b) 


where function / fulfils conditions (A[^-(A|^ (and (A[7]) if needed). 

To close system p.l| ) one needs to define the initial condition 

x{t) = ipiit), y{t) = for t 6 [-T, 0], (flit), ip 2 {t) e C, 


(3.2) 


where C = C([-t, 0],M>) denotes the set of continuous functions defined on the interval [-r,0] 
with non-negative values. Clearly, since in system p.l| ) only y variable has delayed argument model 
dynamics will depend on ^i(O) and ^ 2 ( 0 - 

Theorem 3.1 Let f fulfils assumptions (4[^-(Zl|^, then there exists a global, unique and non-negative 
solution of system ( |3.1| )-( [T2l ). 

Proof : For t e [0, r] y{t - r) = (p 2 {t - t) > 0 is well defined thus for (pi(t), ip 2 (t) e C, under assump¬ 
tion (A[^, there exists an interval [0, f) c [0, r] where the solution (x(t),y(t)) of the system p.l| ) 
exists and is unique, for details see 


Re-writing Eq. p.lb| ) in the integral form (y(t) = ^2(0) exp ^^(1 - x{s))ds) one obtains non¬ 
negativity of y{t) whenever it exists. Now, assume that at point q e (0, f) x(t) starts to be nega¬ 
tive, i.e. x(t) > 0 for 0 < t < ti, x(ti) = 0 and x{ti) < 0. Clearly, for considered f(y) inequality 
i(ti) = d(a) - x(ti) pf{y{t\ - T))x{t\) - fix{t\)y{ti)) = da> > 0, holds. That contradicts the hypothesis 
that solution x(t) starts to be negative for t = ti. 

On the other hand, the non-negativity of x(t) and y{t) imply y(t) < (p2{0)e^\ Moreover, stimulus 
function f(y) fulfils assumptions (A|^-(Aj^. Hence, x < d(a) px) and x{t) < (v^i(O) + “)““■ 

Thus, both variables and theirs derivatives are bounded, so the solution can be extend through the 
point t*. Following the method of steps (see e.g. p^ ) and repeating presented arguments for each 
interval [nr, (n -l- 1)t], n 6 N one completes the proof. ■ 
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3.1 Stability of the steady states 

Clearly, the steady states of the system with delay p.l|) are the same as for the system without de¬ 


lay ( |2.3[ ). However, the system’s stability, already strongly depending on the values of parameters (p 
and 0 ) without delay, now ean also depend on the value of parameter r. 

Theorem 3.2 Let f fulfils conditions (yQ-CQ. The semi-trivial steady state A = (co, 0) of sys¬ 
tem - P.2[ ) is a locally stable node for oj > \ and a saddle for a> < 1 independently of the value 
of parameter r > 0. 

Proof : For r = 0 steady state A is a loeally stable node (for co > 1) or a saddle (for 0 < m < 1) and 
it always exists, for details see p9| . Consider r > 0. Linearising system p.l[ ) around steady state A 
one obtains system without delay whieh is the same as the linearisation of the system p.3| ) without 
delay. Henee, parameter r has no influenee on the loeal stability of semi-trivial steady state A. ■ 

As pointed earlier, for all positive steady states the first eoordinate equals to 1, thus (1,^, ^ > 0 
denotes the positive steady state of the system p.l[). Moreover, to simplify the notation define 


7f ■= Pf(D = p- 


af 


■a-l 


> 0 , 


( 1+^2 

whieh is positive due to the assumption (A|^. 

Theorem 3.3 Let f fulfils conditions (A[^-(A[^, f >Q and jf > 0 is given by p.3[ ). 


(3.3) 


1. If ip - Jf > 0, then the steady state (1, P of system p.l[ ) is unstable for all r > 0. 


2. If\f/-yf < 0, then the steady state (1, of system 63 is locally asymptotically stable for t = 0 
and there exists a critical value t„ such that (1,^) is locally asymptotically stable for t < 

r r 

and is unstable for r > Moreover, for r = T^r, the Hopf bifurcation occurs and periodic 
solutions arise with period equal to In! sjcfi , where 

-(Idpf ■+■ (dcof) ■+■ ^{dcof + Ad^pipcjfi -t- Ald^yff 


= 


— 

‘ cr 




areeos 


¥ + dPf 
d^Jf 


(3.4) 


Proof : Linearising system (3.11 around steady state (1,4") one obtains 


dx(t) 

dt 


= d(-l-t-pfiP) - ip)x(t) - dipyit) pyyyit - r). 


dy(t) 

dt 


-Px(t). 


Note that for steady state (1,4") equality 1 - pf{fi) if/P = co, holds. Thus, the eharaeteristie funetion 
for eaeh positive steady state (1,4”) has the following form 


VF(/l) = -h dcoA - dipp d^yf e' 


-At 


(3.5) 


Clearly, for r = 0 Eq. p.5| ) reads VF(/1) = A^ dcoA - dp{if/ - yf) and for if/ - yf > 0 it has two real 
solutions of opposite signs. Thus, in sueh a ease steady state (1,4") is a saddle. On the other hand, for 
T = 0 and if/-yf < 0 both roots are real and negative or eomplex with negative real parts, henee (1, ^) 
is loeally asymptotieally stable. From the eontinuous dependenee, one knows that for suffieieney 
small delays this result holds, however, the stability of the steady state(s) might ehange for a larger 
values of delay, see e.g. [^|. Sueh a situation is only possible if there exists a eritieal value r = r^r 
and a pair of purely imaginary roots of p.5| ) sueh that roots eross the imaginary axis (in the proper 
direetion) when r exeeeds rir- 
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Let ±W be a pair of pure imaginary roots of ( |3.5| ). Hence, W(W) = -0^+id(o0-d^i//+d^yf(cos6T- 
i sin 6t) = 0 , and 


cos 9t 


m + dif/^ 


> 0, sin Or 


dio6 


> 0 . 


(3.6) 


d^7f d^7f 

A solution of system p. 6 [ ) exists if condition d^{i// - jf) < -6^ is fulfilled. Hence, for 1 /^- 7 / > 0 
the characteristic quasi-polynomial p.5[) of the system p.l|) has no purely imaginary roots and hence 


the real parts of eigenvalues can not change sign for any value of r > 0. That completes the proof of 
part[^ 

Consider ^ - 7 / < 0. Squaring and adding Eqs. p. 6 [ ) one has 

F{e) -.= 6^ + d{2if/^ + du?)e^ + - rfj = O. (3.7) 

Next, substituting rj = 6^ one gets rf + T]({da))^ + Id^iJ/) + (d^)^(d/^ - yj) = 0 with roots given by 

-(2d^if/ + {d(o)^) + ^{2d^if/ + {doSl^y- - Aid^yiifP - y^) 


dm = 


(3.8) 


Clearly, roots p. 8 [) are always real since the expression (2dif/^ + {dooiff- - 4(J4')^(<A^ - 7 ^) is strictly 


positive for if/ - yf < 0. Moreover, for if/ -yf <0 roots have opposite signs. Denote the positive root 


by r]\ > 0 , then there exists 6^ = such that 0 < ( 0 ^)^ -l- 2d^if/ < d^yf, and moreover W(i6^) = 0 
for an infinite sequence of critical values of r given by 


Ty = l 

f 

arccos 

+ 2d^if/' 

+ 2kn 

cr,k 0 ^ 

1 

d^7f 

\ 

/ 


keNo. 


(3.9) 


Denote by r^ the first positive delay fulfilling condition (3.9). To complete the proof one needs to 


show that the real part of eigenvalues of roots cross the imaginary axis at r = with positive speed, 
i.e. ^Re/l(T)|^^^f > 0 and that the re-stabilization of the positive steady state is not possible for larger 


values of r. Here, using the theorem proved by Cooke and Van den Driessche in [^, one can show 
that all the roots of the characteristic quasi-polynomial p.5| ) cross the imaginary axis from the left to 
the right-hand side of the complex plane. The Cooke and Van den Driessche theorem indicates that 
^Re/l(T)|^^^f = hence it is enough to calculate = 4-6^ + 2d9{2if/^ + do/) > 0 for 0 > 0. 

Thus, in the case when the steady state is locally asymptotically stable for t = 0, for r = new 
periodic orbits arise due to the Hopf bifurcation with period equal to 27r/0f. Moreover, destabilized 
steady state (1, p remains unstable for all r > r^r- That completes the proof. ■ 


4 Experimental data, estimation of parameters and model vali¬ 
dation 

The model parameters are obtained by fitting the solutions of the non-scaled version of model p.l| ) 
to the experimental data provided in [ [35| for mice affected by B-cell lymphoma in the spleen (BCLi). 
It is assumed that F is given by the Michaelis-Menten function (F = F 2 defined by ( |2.2| ) with a = 1) 
since it is the most common form of the stimulation function considered in that context in the litera¬ 
ture, see e.g. [|4 21 25 291. Estimation of the ranges for each parameter is based on the experimental 
literature as well as on the assumptions made by other scientists who in particular modelled the BCEi 
and DEBCE (diffuse large B-cell lymphoma) progression or modelled the tumour-immune system 
interactions in general p] T3][T^|^ 2 ^|^ . 
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4.1 Experimental data 

The solutions to non-scaled version of model p.l| ) are fitted to the experimental data for BCLi 
from [ |^ , where the modelled experiment was performed on two sets of animals: the eontrol group 
(BALB/e miee) and ehimerie. Chimerie animals were obtained from BALB/e miee by having, by 
lethal irradiation, their bone marrow destroyed and then having bone marrow transferred from other 
mice to rebuild the immune system. In normal BALB/c mice, one observes that 5x10^ injected BCLi 
cells into the spleens grew rapidly, reaching a plateau after 40 days and causing mice death by 90 
days. In contrast, identically injected groups of chimeric animals demonstrated a different growth 
pattern: the initial rate of growth was similar to the control group, but then it slowed down with tu¬ 
mours reaching their maximum sizes (10 times smaller than in control mice) and a commensurate 
decrease in the number of cells was observed after day 40, see Figure]^ 


4.2 Parameter ranges 


Tumour growth rate (r). Tumour volume duplication time strongly depends on the tumour type. 
For some tumours it can be 20 days, while for others, e.g. colon cancer, it can be years. On the 


other hand, the tumour doubling time for malignant lymphoma is reported to be around 29 days, [37|. 
In [ [3^ , for DLBCL, it was assumed that the tumour doubling time is between 1.4 and 70 days which 
corresponds to a tumour growth rate between 0.0099 and 0.4951 day“' (approximately 0.01 and 0.5 
day“^). That assumption is followed since the interval [0.01,0.5] includes the range 0.05-0.5 proposed 
in |251. Clearly, the tumour growth rate estimated by Kuznetsov et al. in Q (0.18 day"^) also lies in 
that interval. 

Natural lymphocyte death rate (u). The mean lifetime of T lymphocytes from the spleen and 


the blood was estimated to be approximately 30 days or more, [^. Thus, the degradation rate, which 
is l/(mean life time), is assumed to be in the range 1.25 x 10“^ - 5 x 10“^ day“\ which corresponds to 
20 - 80 days mean life time. The values of the death rate 0.03125 and 0.0412 day“^ estimated in |251 
and Q, respectively, also lie in the assumed range. 

Immune response function (F, a, k 2 , a). It is assumed that the immune system response term 
has the same form as those considered e.g. in [|4}[T3||2Tj|25j|29| or [[3^, that is F = Y°'+ Y") 
with a = 1. In Q the value of the steepens of the immune response function (^ 2 ) was estimated to be 
2.019 X 10^ cells, while in [251 the value 1.5109 x 10^ cells. In this study the range of 10’ - 2 x 10^ 
cells for ^2 is chosen. To determine the range of the maximal immune system response rate (a) [ |25| 
is followed and range 0-2.5 day“', which contains the value 0.1245 day“^ estimated in Q, is taken. 

Tumour independent production rate of effector cells (w). Following [|25| the range for the 


constant influx of effector cells (w) from the immune system in the absence of tumour is chosen to 
be 3.2 X 10^ - 3.2 x 10"* cells day“k The value 1.3 x 10"^ cells day“' estimated by Q belongs to the 
considered range. 

Tumour-induced inactivation rate of effector cells (b). In [ |25| the range for the fraction of 
immune cells inactivated due to the interactions with tumour cells (w) was assumed to be 10“'^ - 5 x 
10“’ cells"' day"'. The value 3.422 x 10"'" cells"' day"' estimated in Q belongs to the assumed 
range. 

Effector induced elimination rate of tumour cells (c). In [ [T8| the range for the fraction of cancer 
cells killed per effector cell was assumed to be 10"" - 10"’ cells"'day"'. However, in Q the value of 
that parameter was estimated to be equal to 1.101 x 10"’ cells"'day"'. Hence, the range is extended 
and it is assumed to be c 6 [10"", 3 x 10"’]. 

Time delay in the immune system activation (t). In this study it is assumed that time delay in 
immune system response is in range from Ih till 5 days, i.e. 0.0416-5 days. 

All ranges of the model parameters together with corresponding references are given in Table 
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Table 2: Parameter ranges with refereneed used for the fitting proeedure to the experimental data 
from [[^ ■ _ 


Par. 

Biological interpretation 

Unit 

Range 

Reference 

w 

tumour independent production 
rate of effector cells 

cells day“^ 

3.2x 10^ -3.2x 10^ 

\25\,^ 


u 

natural lymphocyte death rate 

day“^ 

1.25 X 10-2-5 X 10-2 

U) 

00 

1 K> 

51, [4|. 

a 

maximal immune response rate 

day“^ 

0-2.5 

[251, ^ 

H 

b 

fraction of immune cells inacti¬ 
vated in interactions with tumour 
cells 

cells" ^ day"^ 

10-12 - 5 X 10-2 

[251, [^ 

g 

r 

tumour growth rate 

day"^ 

0.01-0.5 

[361, [ 
[37! 

25uir 

c 

fraction of cancer cells killed per 
effector cell 

cell"^ day"^ 

10-11 - 3 X 10-2 

[18L| 


f 

time delay 

day 

0.0416 - 5 

ass. 


ki 

steepens of immune response 

cells 

102 - 2 X 10*1 

[251, [^ 

l\ 

a 

detrmines type of immune re¬ 
sponse 

— 

1 

4|, [1 
[21|,[2 

31, [291, 
5|,[36[ 


4.3 Estimation of model parameters 

Clearly, if there are no tumour eells in the system (7 = 0), then immune-activity stays at the constant 
level calculated directly from ( |2.1a[ ). Hence, the initial condition for the variable X is a constant 
equal to w/u. On the other hand, to determine the level of the tumour cells at the beginning of the 
experiment (at s = 0) one can directly use the data from [35| taking 02(0) = 5 x 10^ cells. Thus, the 
initial condition for non-scaled version of p.ll) reads 


w 

Oi(5) = — for 5 6 [-T,0], and 02 ( 5 ') 
u 


|o for 5 6 [-T, 0), 

js X 10^ for 5 = 0, 


(4.1) 


and due to the dependencies on the parameters w, u and r = r/r the initial condition is estimated 
within the fitting procedure. 

The fitting procedure is performed in the commercial MATLAB® computational environment 
using the PSO function. PSO finds a minimum of a function of several variables using the particle 
swarm optimization algorithm based on the stochastic optimization technique that was inspired by 
the social behaviour of bird flocking or fish schooling. The algorithm was originally introduced in 
1995 by Kennedy and Eberhart, [ |^ , and later extended by Shi and Eberhart, [ |40| . Nowadays it has 
a number of extensions, however the main idea is the same a swarm of particles moves around in 
the search space and the movements of the individual particles are influenced by the improvements 
discovered by the others in the swarm. In [411 Clerc and Kennedy introduced a constriction factor 
in PSO, which was later shown to be superior to the inertia factors. That version of algorithm is 
used here. Solutions of delay differential equations (DDEs) are computed using dde23 function with 
tolerance RelTol =le-7, and AbsTol=le-l 1. 

In the presented study the objective function is the standard one based on the calculation of mean 
square error (MSE): 


MSE 


.. m 

-y 

m 

(=1 


err- 


erri 


^ata _ 

^data ’ 


where m is the number of experimental measurements, are the measured values and are the 
values of the solution of DDEs for respective time points. 
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Tables 

system 


Fitted non-scaled model parameters, steady states and theirs stability for r = 0 for non-scaled 


3.11. Parameter a= 1. 



control 

chimeric 

parameters 

w 

3.0337 X 10^ 

u 

0.0196 

a 

0.1452 

0.2835 

ki 

1.3604 X 10^ 

1.8720 X 10^ 

b 

2.1010 X 10-^° 

6.6828 X 10-‘° 

r 

0.4695 

c 

1.7608 X 10-’ 

2.0761 X 10-’ 

f 

3.1869 

0.0581 

steady states and theirs stability for t = 0 

semi trivial unstable 

(1.5472 X 10^0) 

(1.5472 X 10^0) 

positive stable 

(2.6666 X lOM.0548 X 10^) 

(2.2616 X 10^7.8919x 10^) 

positive unstable 

(2.6666 X 10^ 5.0530 X 10^) 

(2.2616 X 10^2.1986x 10^) 


Table 4: Fitted model parameters for system (3.1), parameter q'=1 



d 

to 

P 


T 

control 

0.0418 

0.5802 

7.4046 

1.4576 

1.4964 

chimeric 

0.6841 

14.4583 

6.3802 

0.0273 


First, the solutions of the non-scaled version of system p.l| ) are fitted to the data from the chimera 
experiment for BCLj tumour cell population in the spleen of the BALB/c mice over 110 days (de¬ 
scribed in [^) with the initial condition given by ( |4.1[ ) that depends on w, u and t parameter values. 
The result of the fitting procedure is presented in Figure (left) showing a good agreement with the 
experimental data (with MSB at the level of 3.57%) and thus positively validates the proposed model. 

Next, the parameters w, u and r are fixed (as for chimera experiment) and the rest of the parameters 
are fitted to the 90 day control mice experiment to further validate the model. This time a very good 
agreement with the experimental data with MSE at the level of 0.14% is obtained, see Figure]^ (right). 
The corresponding values of the estimated parameters for both experiments are given in Table 


5 Results 

Comparing the sets of parameters (see Table for both experiments estimated via the fitting proce¬ 
dures one immediately sees that for the chimeric set of parameters parameter a is two times larger 


Table 5: Summary of the bifurcation phenomena for scaled parameters given in Table and r treated 
as a bifurcation parameter, a = 1. 


Param. type 

steady state 

stability r = 0 

Ar 

period 

bif. type 

control 

(0.5802, 0) 

unstable 

unstable for r > 0 

(1,0.0775) 

stable 

1.1776 

50.1363 

supercritical 

(1,3.7145) 

unstable 

unstable for r > 0 

chimera 

(0.6841, 0) 

unstable 

unstable for r > 0 

(1,0.0422) 

stable 

1.2227 

57.3715 

supercritical 

(1, 1.1745) 

unstable 

unstable for r > 0 
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Figure 2: (Left) Result of the fitting proeedure of the trajeetories of the non-sealed version of sys¬ 
tem p.l[ ) to the experimental data for BCLj tumour eell population in the spleen of the ehimerie miee 
over 110 days experiment deseribed in p5| . MSE=0.0357. (Right) Result of the fitting proeedure of 
the trajectories of the non-sealed version of system p.l[ ) to the experimental data for BCLi tumour 
cell population in the spleen of the BALB/c mice over 90 days control experiment described in [ f35| . 
MSE=0.0014. In both panels solutions of the non-sealed version of system p.l| ) are computed with 
RelTol=le-7 and AbsTol=le-ll. All used parameter values, for both panels, are given in Table 
Bars indicate the measurement errors as reported in [[35|. 


than for the control set. Thus, together with slightly larger k 2 (for the chimeric experiment) implies 
a stronger stimulation of the immune system due to the presence of the tumour cells. Additionally, 
the time lag of that response is much shorter, approx. 1.4 h, for the chimeric set of parameters than 
for the control, approx. 76.5 h, indicating a more effective immune system reaction in the presence 
of the tumour cells for the chimeric mice. Eor both the chimeric and control sets, the estimated value 
of the tumour cells killed per effector cells (c) is very similar. Despite the large range for parameter b 
(10“^^ - 5 X 10“^, see Tableand that no further assumptions on the range of b to differentiate the 
chimeric and control experiments are made, the fitted values obtain the same order of magnitude. 

The simulations presented in Eigure again validate positively the model and confirm results 
from [351 showing that for larger than control initial injection of populations of the tumour BCEi 
cells into the spleens of the chimeric mice leads to the same growth pattern. Experiments performed 
in [351 indicate that even if one injects 5 x 10^ cells after around 110 days, due to the immune system 
response, the number of tumour cells is reduced to the almost the same amount as for 5x10^ initially 
injected cells. To reproduce that phenomena all used parameter values were fixed, as given in Table 
(column chimeric) and only a part of the initial condition that is 02(0) was changed, see Eigurej^Erom 
the mathematical point of view such behaviour shows that the stable steady state, that can be identified 
as a dormant tumour state, has large enough basin of attraction. 

Eor the sets of estimated parameters that have been scaled and are given in Table|^the system ( [3.1[ ) 
has one semi-trivial steady state which is equal to (0.5802,0) and (0.6841,0) for the control and 
chimeric parameter sets, respectively. Eor both sets of parameters o) <\, thus both semi-trivial steady 
states are unstable. Moreover, according to Theorem [3^ they stay unstable independently of the value 
of time delay. 

Additionally, for both sets of parameters, assumptions of Proposition [2.2[ are fulfilled and thus, 
for both cases the system p.l[ ) has two positive steady states: = (1,0.0775) and Cc = (1,3.7145) 

for control; Bchim = (1,0.0422) and Cchim = (1,1.1745) for chimera. Interesting is that the values of 
the positive steady states Bchim and Be are similar. Eor r = 0 one can easily check that the steady 
states Be and Bchim are locally asymptotically stable, while the steady states Cc and Cchim are unstable, 
see Table 1^ Considering r as a bifurcation parameter and using results of Theorem 3.3 one obtains 
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Figure 3: Trajectories of the non-scaled version of tumour-immune system ( |3.1| ) for initial condition 
as defined in ( |4.1| ) with different values of 02(0): 5 x 10^ (dash green line), 5 x 10® (dash-dot blue 
line) and 5 x 10^ (solid red line). All used parameter values are given in Table (column chimeric), 
5 6 [0,110]. All solutions computed with RelTol=le-7 and AbsTol=le-ll. 




Figure 4: (Left) Critical values of the bifurcation parameter as the functions of to (left) and as the 
functions of p (right) for both chimeric (red dash line) and control (blue solid line) sets of parameters. 
All used parameter values (except oj, p and r) are given in Table The stability regions lie below the 
solid blue and the dashed red lines. The horizontal grey dash-dot line in the left panel indicates the 
value of delay estimated via the fitting procedure for the control experiment r = 1.4964. 


that independently of the value of time delay the steady states and Cchim remain unstable. The 
situation differs for steady states Be and Behwu whose stability changes for sufficiently large values 
of parameter r. Moreover, Theorem |3.3| indicates that there can be only a single stability switch for 
each Be and Behim steady states for time delay treated as a bifurcation parameter. For the considered 
sets of parameters the critical values of the bifurcation parameter are 1,1776 and 1.2227 for control 
and chimeric parameter sets, respectively. Theorem |3.3| and the results of the numerical analysis 
indicate that for sufficiently large time delay, the periodic behaviour of the solutions of system p.lj ) 
should be observed if the initial condition is properly chosen.Clearly, for the estimated chimeric set 
of parameters the delay value dose not exceed the critical value, thus for that set of parameters the 
steady state Behim is locally asymptotically stable. For the control set of parameters. Be is unstable 
since for r = 1.1776 the Hopf bifurcation takes place. 

The functional dependences of the critical values of the bifurcation parameter {Ter) for the steady 
states Be and Behim on the model parameters oj and p are presented in Figure Q where all parameters 
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are fixed (as in Table except oj, p and r. Figure left indicates that the critical values of time 
delays, for both chimeric and control sets of parameters, are the increasing functions background 
immunity. Thus, the increase of the native immunity increases the stability region of the steady states 
Be and Behim- Proposition |2.2| indicates that for 0 < m < 1 (that is the case of the estimated values of 
parameters) and p ^ (+ Vl - (that is p ^ 9.5355 and p ^ 3.4419, for chimeric and control 
sets of parameters, respectively) the positive steady states merge, thus the critical values can only 
be computed for sufficiently large values of p. Additionally, there are no upper limits for p. Similarly 
as in the previous case. Ter for both sets of parameters are increasing functions of p. However, this 
time they are concave functions. Thus, increasing the stimulation strength enlarges the stability region 
of the steady states Be and Behim, however in a different manner. 

In the general case a large number of model parameters makes the analytical study of the nature 
of the Hopf bifurcations, i.e. of the stability of the arising periodic solutions, hard to perform. Thus, 
for the estimated model parameters the BIFTOOL package (a Matlab package for numerical bifurca¬ 
tion analysis for DDEs developed by Engelborghs et al.) is used, for details and the documentation 
see p2| , p3| or visit http://twr.cs.kuleuven.be/research/software/delay/ddebiftool.shtml, InEigure]^ 
the amplitudes of the periodic solutions as a function of the bifurcation parameter r for the control 
and chimeric sets of parameters (see Tablefor model ( |3.1[ ) are presented. Clearly, from the figure 
one reads that periodic solutions are born for critical values of r given in the Table and exist for t 
greater than the critical values of the bifurcation parameter, indicating that the Hopf bifurcations are 
supercritical. In both cases, the amplitudes are comparable, having the same order of magnitude. This 
result holds for the non-scaled model as well since the parameter r is held constant across both cases, 
whilst the fitted value of c has the same order of magnitude between cases. 



Figure 5: Amplitudes of the periodic solutions arising due to the Hopf Bifurcation as a function of 
the bifurcation parameter r for chimeric (left) and control (right) sets of parameters given in Table 
for model (|3.1|). In green the predicted values are marked, while in blue the corrected ones. 


In Eigure|^(top) the period of periodic solutions arising from the Hopf bifurcations as a function 
of T for control and chimeric sets of parameters are presented. The directly calculated (Theorem 3.3) 
values of the periods for critical values of bifurcation parameter are 50.1363 and 57.3715 for control 
and chimeric sets of parameters, respectively, see Table Calculating the same for the non-scaled 
model gives periods of approx. 107 and 122 days for the control and chimeric sets of parameters, 
respectively. These pattern agrees with 3-4 months oscillations observed for the clinical data for 
certain kinds of human leukemias and also agrees with oscillatory growth of human and rat tumour 


cell lines reported in [44-461 or in [471 for non-Hodgkin’s lymphoma. 


Interesting is that such recurrent patterns of tumour remission and re-growth have also been clin¬ 
ically observed for immunized BAEB/c mice, where mice where initially immunized with BCEi- 
Idiotype before the injection of 10^ BCEi tumour cells and without any further immunization, [35 j. 
Eor such experiment the level of Anti-Id in serum and spleen sizes where followed for a period of 115 
days showing that in several mice the spontanous regression followed by splenomegaly was observed. 
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Figure 6: Top: Periods of periodic solutions arising from the Hopf bifurcations as a function of 
bifurcation parameter r for chimeric set of parameters (left) and control set of parameters (right) for 
model p.l| ). Middle and bottom: Floquet multipliers plotted on the complex plain for the periodic 
solutions near the bifurcation points for chimeric (left) and control (right) sets of parameters for 
model p.l[ ). In the bottom row the zoom of the results form the middle row is presented. For all 
simulations model parameters are given in Table 


compare Figure 5 in [ [35| |. Moreover, Figure]^ (top) indicates that the periods of the periodic solutions 
increase with the increase of the bifurcation parameter thus, for larger time lags in the reaction of 
the immune system to the presence of antigens, larger periods of oscillation of the tumour mass and 
immunity are observed. 


In Figure (middle and bottom) the Floquet multipliers plotted on the complex plain for both 
sets of estimated parameters are presented. Clearly, the multiplier (0,1) is the trivial one that always 
exists for the autonomous systems, while the ones inside the circle has modulus smaller than 1. Since 
periodic solutions are stable if all, except (0,1), multipliers have modulus smaller than one, one con¬ 
cludes that for the parameter values given in Table the arising periodic solutions are stable. This 
implies that the supercritical Hopf bifurcations for both sets of parameters are observed, confirming 
our previous results presented in Figure]^ From the biological or medical point of view supercritical 
Hopf bifurcation is considered as a safe bifurcation (of course assuming that the amplitude of periodic 
solutions can not be too large), since although the steady state becomes unstable for larger values of 
the bifurcation parameter the new born periodic solutions are stable. In such situations a local stable 
limit cycle around the dormant tumour steady state occurs and again some kind of boundedness of tu¬ 
mour growth is observed. Such oscillations around the steady state corresponds to situation where the 
tumour grows and shrinks with no application of additional treatment, see the experimental tumour 
volume data from [[44l|45l or p)6|. 
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6 Discussion 


In this paper the interaetion between tumour and immune system by a system of non-linear differential 
equations with diserete delay, representing time lag in the reaetion of immune system to reeognition 
of the antigens, is described. The mathematical properties of proposed model such as: existence, 
uniqueness and non-negativity of solutions are studied. The presented study focusses on the on the 
stability of the steady states depending on the value of the bifurcation parameter, which is the time 
delay. The direct conditions for the stability of the steady states are derived in this paper. In the 
general case, it is shown that depending on assumptions on the model parameters the stability of the 
system changes due to Hopf bifurcation for sufficiently large delay. Moreover, only one such stability 
switch is possible. Additionally, it is shown that the unstable (for r = 0) steady states remain unstable 
independently of the values of delay. 

Next, the proposed model is calibrated with two sets of the experimental data for BCLi tumour 


cells injected to the mice spleens, [351, with MSEs at the level of 3.57% and 0.14%. For the estimated 
values of parameters, the proposed model is further validated with other experiments for different 
initial numbers of injected cells. For both sets of estimated values of the parameters the considered 
model has two positive steady states and one unstable semi-trivial steady state. For the chimeric set 
of parameters there are stable and unstable positive steady states, which indicate a local stabilization 
(dormancy) of tumour growth at finite size if the initial size of tumour is sufficiently small (as in 
case of considered experimental setup) or infinite tumour growth otherwise. For the estimated control 
set of parameters both positive steady states are unstable, however the periodic orbits exist. That 
means that depending on the initial conditions the unrestricted tumour growth (as for the control set 
of parameters) or periodic oscillations of tumour size are observed. 

The second coordinates of the positive stable steady state for the chimeric set of the parameters 
and positive steady state (that is stable for sufficiently small delay) for the control set of parameters 
are similar. Hence, by decreasing the time response of the immune system to the presence of the 
tumour cells one might stabilize the uncontrolled tumour growth in case of the control experiment. 

Additionally, for both estimated sets of parameters (except time delay treated as a bifurcation 
parameter) the presented analytical study is extended by investigating the stability of the new born 
periodic orbits and thus the type of existing Hopf bifurcations. It should be pointed out that the 
oscillatory behaviour of the multicellular tumour volume was experimentally observed, even with¬ 
out administration of any treatment, as reported for the in vitro experiments for human and rat cell 


lines reported in and [441. Such oscillations are explained e.g. as the effect of non-linear feed¬ 
back between proliferation at the multicellular tumours surface and invasion of the surrounding by 
tumour cells leaving the spheroid surface however other scenarios are also discussed. The oscillating 
growth phenomenon is also observed for the in vivo studies and take place for e.g. the mice initially 
immunized with BCFi-Idiotype before the injection of 10^ BCFi tumour cells or in non-Hodgkin’s 


lymphoma [47[. 


Although the model is rather simple, its dynamics might be complex and reflect a number of bio¬ 
medical phenomena such as: therapy failure considered as a modification of model parameters, i.e. 
an unrestricted growth of the tumour; stabilisation (at least temporary) of tumour growth at a finite 
size - dormancy; and desirable regression of the tumour. If the initial size of the tumour (or rather 
amount of antigens) is strictly greater than zero, then according to the values of the model parameters 
seven different scenarios are possible: (i) asymptotic unrestricted growth of the tumour independently 
of the initial tumour size; (ii) stabilisation of the tumour growth independently of initial tumour size 
(existence of globally stable equilibrium); (iii) global eradication of the tumour by the immune system 
independently of the initial tumour size; (iv) local stabilisation of tumour growth at a certain finite 
size (if the initial tumour size is sufficiently small one observes stabilisation of the tumour growth 
at a finite size, otherwise unrestricted growth); (v) local eradication of the tumour by the immune 
system (for sufficiently small initial size tumour will be eradicated by the immune system, otherwise 
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it will grow unrestrictedly); (vi) bi-stability with eradication for sufficiently small initial tumour, 
stabilisation on a finite size of tumour of medium initial size and unrestricted growth for large initial 
tumour when two locally stable steady states co-exist; and finally (vii) recurrent patterns of tumour 
remission and re-growth. The first six of these features can be observed for the model without time 
delay, while the last one it is not possible since the ODE model has no periodic solutions. Hence, 
the DDE model can reflect the oscillatory phenomenon observed for a range of tumour types and in 
particular in lymphoma. 

The main goal of the any-anticancer therapies is to treat disease or at least control its progression. 
Eorm the mathematical point of view this means that one wishes to ensure that the trajectories of 
the system are in the region of attraction of the proper steady states. To do that at the level of the 
considered model one can try to increase or decrease some of the parameters. In the context of 
immunotherapy several attempts to increase the immune system efficiency or the creation of anti¬ 
tumour immunity were already proposed in many different ways [[T|[^^^. One can consider the 
increase of the native immunity e.g. by proper incubation of a patient’s natural killer lymphocytes to 
enhance their activity and later by their infusion of them into the patient body. In the mathematical 
language that would mean the increase of the co parameter which is crucial for reaching the global 
or local eradication of the tumour or the bi-stability phenomenon (depending on the values of the 
other system parameters). Recall that a> = wcl(ur). Clearly, since the natural lymphocyte death rate 
(u) and the tumour growth rate (r) are rather out of reach of immunotherapy one should thus focus 
on increasing w and c parameters, that is stimulation of the production strength and the effectiveness 
of effector cells. The increase of the effective cells can be reached by the infusions of the cells 
grown ex vivo, while the increase of the killing effectiveness of the cells can be obtained as the 
result of the infusion of properly chosen cytokines, O' On the other hand, parameter r can be for 
example decreased by different types of chemotherapies. Another parameter that has an essential role 
in determining the range of the controlled tumour growth is the stimulation strength p = a/u. If 
the tumour antigens are not presented to the lymphocytes, one can increase a by method of artificial 
loading of the of antigen-presenting calls with tumour antigen. The effectiveness of such an approach 
will be dictated by the values of the background immunity and strength of the anti-immune activity. 
In this paper, it is shown that for the model with discrete delay and estimated parameters an increase 
of the background immunity (m), as well as the stimulation strength (p) of the tumour cells on the 
immune system, causes the increase of the critical delay value for which the Hopf bifurcation is 
observed, enlarging the stability region for the positive dormant steady state, see Eigurej^ Hence, by 
increasing the background immunity (range is reachable) one might stabilize the unrestricted tumour 
growth observed for the control experiment. 

On the other hand, if one wishes to model the effectiveness of a given therapy or try to increase a 
given therapeutic effect one should rather directly include these elements in the system, e.g. by adding 
new terms and considering the optimization problem. That of course would force the additional 
calibration of the extended model with additional sets of experimental data required. Researchers who 
are interested in collaborating with experimentalists or clinicians who have access to patients or in 
vivo animal trial data should feel encouraged to follow this idea. Recently the very promising concept 
of an anti-tumour vaccines has arisen, nura- However, due to the discrete nature of the vaccination 
program, the presented model would need to be adapted to include an impulse-response framework. 
In that context, the presented calibrated model can be a base for the future study of the effectiveness of 
a different types of therapies, including the chemotherapy or combined chemo-immunotherapy. One 
should feel free to modify the proposed model by fitting it to appearing needs. Eor example one can 
consider instead of the exponential tumour growth rather the Gompertz one, which takes into account 
the finite environmental capacity for tumour cells. That might cause an essential qualitative difference 
between models. However, it is still an open question. On the other hand, to fit even better proposed 
model with the experimental data one could instead of the discrete delay consider the distributed one 
that often better reflects the real measurements. 
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